1 Wymagania

1.2 Executive summary

raport powinien zaczynać się od rozdziału podsumowującego całą analizę, streszczającego najważniejsze spostrzeżenia analityka

1.3 Lista wymagan minimalnych [kolejnosc dowolna]

Kod wyliczający wykorzystane biblioteki. Kod zapewniający powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych. Kod pozwalający wczytać dane z pliku. Kod czyszczący dane (np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych). Sekcję podsumowującą rozmiar zbioru i podstawowe statystyki. Szczegółową analizę wartości atrybutów. Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji. Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie. Sekcję próbującą stworzyć klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji. Analizę ważności atrybutów najlepszego znalezionego modelu. Dodatkowo punktowane będzie wykonanie analizy typowej dla danych klinicznych, np. regresji logistycznej wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors) lub regresji Coxa (ang. Cox Proportional-Hazards Model).

1.4 Dodatkowe uwagi

Analityk nie musi, a nawet nie powinien, ograniczać się do powyższych punktów. Wszelkie dodatkowe techniki analizy danych, wizualizacje, spostrzeżenia będą pozytywnie wpływały na ocenę.

Ewentualne konkluzje, znalezione zależności warto potwierdzić dokonując sprawdzenia istniejących wyników badań w literaturze naukowej (np. na Google Scholar czy PubMed).

2 Wykonanie

2.1 Podsumowanie

TBA

2.2 Import bibliotek

library(xlsx)
library(DT)
library(knitr)
library(dplyr)
library(tidyr)
library(janitor)
library(imputeTS)
library(lares)
library(plotly)
library(caret)

2.3 Wczytanie i wstepna analiza danych

2.3.1 Wczytanie danych,

filename <- 'res/wuhan_blood_sample_data_Jan_Feb_2020.xlsx'
raw_data <- read.xlsx(filename, 1)
raw_data <- as_tibble(raw_data)

2.3.2 Sprawdzenie surowych danych

# OK
#how many NA in variables
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596
# var min i max
min(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] -1
max(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)
## [1] 50000
glimpse(raw_data)
## Rows: 6,120
## Columns: 81
## $ PATIENT_ID                                                    <dbl> 1, NA...
## $ RE_DATE                                                       <dttm> 2020...
## $ age                                                           <dbl> 73, 7...
## $ gender                                                        <dbl> 1, 1,...
## $ Admission.time                                                <dttm> 2020...
## $ Discharge.time                                                <dttm> 2020...
## $ outcome                                                       <dbl> 0, 0,...
## $ Hypersensitive.cardiac.troponinI                              <dbl> NA, N...
## $ hemoglobin                                                    <dbl> NA, 1...
## $ Serum.chloride                                                <dbl> NA, N...
## $ Prothrombin.time                                              <dbl> NA, N...
## $ procalcitonin                                                 <dbl> NA, N...
## $ eosinophils...                                                <dbl> NA, 0...
## $ Interleukin.2.receptor                                        <dbl> NA, N...
## $ Alkaline.phosphatase                                          <dbl> NA, N...
## $ albumin                                                       <dbl> NA, N...
## $ basophil...                                                   <dbl> NA, 0...
## $ Interleukin.10                                                <dbl> NA, N...
## $ Total.bilirubin                                               <dbl> NA, N...
## $ Platelet.count                                                <dbl> NA, 1...
## $ monocytes...                                                  <dbl> NA, 1...
## $ antithrombin                                                  <dbl> NA, N...
## $ Interleukin.8                                                 <dbl> NA, N...
## $ indirect.bilirubin                                            <dbl> NA, N...
## $ Red.blood.cell.distribution.width.                            <dbl> NA, 1...
## $ neutrophils...                                                <dbl> NA, 6...
## $ total.protein                                                 <dbl> NA, N...
## $ Quantification.of.Treponema.pallidum.antibodies               <dbl> NA, N...
## $ Prothrombin.activity                                          <dbl> NA, N...
## $ HBsAg                                                         <dbl> NA, N...
## $ mean.corpuscular.volume                                       <dbl> NA, 9...
## $ hematocrit                                                    <dbl> NA, 3...
## $ White.blood.cell.count                                        <dbl> NA, 3...
## $ Tumor.necrosis.factorÎ.                                       <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.concentration                     <dbl> NA, 3...
## $ fibrinogen                                                    <dbl> NA, N...
## $ Interleukin.1β                                               <dbl> NA, N...
## $ Urea                                                          <dbl> NA, N...
## $ lymphocyte.count                                              <dbl> NA, 0...
## $ PH.value                                                      <dbl> 7.415...
## $ Red.blood.cell.count                                          <dbl> NA, 4...
## $ Eosinophil.count                                              <dbl> NA, 0...
## $ Corrected.calcium                                             <dbl> NA, N...
## $ Serum.potassium                                               <dbl> NA, N...
## $ glucose                                                       <dbl> NA, N...
## $ neutrophils.count                                             <dbl> NA, 2...
## $ Direct.bilirubin                                              <dbl> NA, N...
## $ Mean.platelet.volume                                          <dbl> NA, 1...
## $ ferritin                                                      <dbl> NA, N...
## $ RBC.distribution.width.SD                                     <dbl> NA, 4...
## $ Thrombin.time                                                 <dbl> NA, N...
## $ X...lymphocyte                                                <dbl> NA, 2...
## $ HCV.antibody.quantification                                   <dbl> NA, N...
## $ D.D.dimer                                                     <dbl> NA, N...
## $ Total.cholesterol                                             <dbl> NA, N...
## $ aspartate.aminotransferase                                    <dbl> NA, N...
## $ Uric.acid                                                     <dbl> NA, N...
## $ HCO3.                                                         <dbl> NA, N...
## $ calcium                                                       <dbl> NA, N...
## $ Amino.terminal.brain.natriuretic.peptide.precursor.NT.proBNP. <dbl> NA, N...
## $ Lactate.dehydrogenase                                         <dbl> NA, N...
## $ platelet.large.cell.ratio.                                    <dbl> NA, 3...
## $ Interleukin.6                                                 <dbl> NA, N...
## $ Fibrin.degradation.products                                   <dbl> NA, N...
## $ monocytes.count                                               <dbl> NA, 0...
## $ PLT.distribution.width                                        <dbl> NA, 1...
## $ globulin                                                      <dbl> NA, N...
## $ γ.glutamyl.transpeptidase                                    <dbl> NA, N...
## $ International.standard.ratio                                  <dbl> NA, N...
## $ basophil.count...                                             <dbl> NA, 0...
## $ X2019.nCoV.nucleic.acid.detection                             <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.                                  <dbl> NA, 3...
## $ Activation.of.partial.thromboplastin.time                     <dbl> NA, N...
## $ High.sensitivity.C.reactive.protein                           <dbl> NA, N...
## $ HIV.antibody.quantification                                   <dbl> NA, N...
## $ serum.sodium                                                  <dbl> NA, N...
## $ thrombocytocrit                                               <dbl> NA, 0...
## $ ESR                                                           <dbl> NA, N...
## $ glutamic.pyruvic.transaminase                                 <dbl> NA, N...
## $ eGFR                                                          <dbl> NA, N...
## $ creatinine                                                    <dbl> NA, N...
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770596

2.4 Transformacja danych

(np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych)

# NOT OK
#ogarnac te brakujace daty -- tyle ze ich nie ma?
# kolumny z samimy na


#replace -1 with NA
raw_data[raw_data==-1]<-NA

#PATIENT_ID
id_filled <- raw_data %>% fill(PATIENT_ID)


#OTHER NA VALUES
#is there any row that doesn't have NA value?
full_row_count <-dim(id_filled[complete.cases(id_filled),])[1]
print(full_row_count)
## [1] 0
# NA values in variables:
mean(is.na(id_filled %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8770748
# there are two patients with negative days, as their only blood sample results arrived one day after their clinical outcome


#remove rows where all variables are empty
vars <- colnames(id_filled)[-(1:7)]
no_empty_rows<- id_filled[rowSums(is.na(id_filled[vars])) != length(vars), ]
no_empty_cols <- no_empty_rows[colSums(!is.na(no_empty_rows)) > 0]
# NA values in variables:
mean(is.na(no_empty_rows %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))
## [1] 0.8658041
colnames_cleaned <- no_empty_cols %>% clean_names()

NA Values:

clean_NA<-function(column){
  not_NA_count<-sum(!is.na(column))
  if (not_NA_count>=2){
    column <- na_interpolation(column, option = "linear")
    column
  }

  else if (not_NA_count==1){
    val <- first(na.omit(column))
    column[is.na(column)] <- val
    column
  }
  column
}

cleaned<- colnames_cleaned%>% group_by(patient_id) %>% mutate_each(funs(clean_NA))
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `mutate_each_()` is deprecated as of dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
var_df<-cleaned[-(1:7)]

2.5 Prezentacja czystych danych

“Sekcja podsumowująca rozmiar zbioru i podstawowe statystyki” Szczegółowa analiza wartości atrybutów.

DT::datatable(cleaned, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))
## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html
all_dim <-dim(cleaned)
tmp <- cleaned %>% select(patient_id, outcome, gender) %>% group_by(patient_id) %>% summarise(outcome_count = mean(outcome), gender_count =  mean(gender))
## `summarise()` ungrouping output (override with `.groups` argument)
# ile pacjentow
patients_count <- length(distinct(cleaned, patient_id)$patient_id)
# ile pomiarow
measurements_count <- all_dim[1]
# ile srednio pomiarow na pacjenta
mean_measures <- measurements_count/patients_count
# ile kobiet/mezczyzn
genders <- tmp %>% count(gender_count) #Male 224, Female 151 # WYKRES
# ile umarlo/przezylo
outcomes <- tmp %>% count(outcome_count) #201 recovered, 174 died # WYKRES
# ile kolumn
columns_count <- all_dim[2]
# ile atrybutow
vars_count <- all_dim[2]-7
# ile zostalo wartosci NA
na_left <- round(100*mean(is.na(cleaned)), 0)

titles <- c("Liczba pacjentów", "Liczba pomiarów", "Średnia liczba pomiarów", "K", "M", "Smierc", "Zycie :(", "Liczba wierszy", "Liczba zmiennych", "Brakujace wartosci [%]")
values <- c(patients_count,measurements_count,mean_measures,1, 2, 0,1, columns_count,vars_count, na_left)
info_table <- data_frame(
  titles,
  values
)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
knitr::kable(info_table)
titles values
Liczba pacjentów 361.00000
Liczba pomiarów 5606.00000
Srednia liczba pomiarów 15.52909
K 1.00000
M 2.00000
Smierc 0.00000
Zycie :( 1.00000
Liczba wierszy 80.00000
Liczba zmiennych 73.00000
Brakujace wartosci [%] 6.00000 
cleaned %>% distr(outcome, gender)

#wykres death_bothgenders, alive_bothgenders
invisible(remove(tmp))

plot_data <- cleaned%>%select(patient_id, gender, admission_time, discharge_time, outcome) %>%group_by(patient_id) %>% summarise_all(funs(first))
admission_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    facet_wrap(~outcome) +
    layer(data=plot_data, 
          mapping=aes(
              x=admission_time, 
              y=discharge_time, colour=factor(gender)
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(admission_plot)

2.6 Analiza wartosci atrybutow

Szczegółową analizę wartości atrybutów.

#data.frame(unclass(summary(tmp)), check.names = FALSE, stringsAsFactors = FALSE)
res <- var_df%>% sapply(function(x){c(min=min(x,na.rm = TRUE), median=median(x,na.rm = TRUE), mean = mean(x,na.rm = TRUE),max=max(x,na.rm = TRUE), variance=var(x,na.rm = TRUE), standard_deviation=sd(x,na.rm = TRUE), na_count = sum(is.na(x)))}) %>%  data.frame()
rd<- dim(res)
rows<-rd[1]
cols<-rd[2]
plot_data<-data.frame(c(1:cols),c(1:cols))

DT::datatable(res, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))
#plt<-ggplot(gather(var_df), aes(value)) + 
#    geom_histogram(bins = 10) + 
#    facet_wrap_paginate(~key,nrow = 1, 
#                     ncol = 1, 
#                     scales = "free")
#ggplotly(plt)

page_plot<-function(page_no, prow=3, pcol=3){
  print(
  ggplot(gather(var_df), na.action="na.omit", aes(value)) +   ggtitle(" title")+
    geom_histogram(bins = 10) + 
    facet_wrap_paginate
  (~key, ncol = pcol, nrow = prow, scales = "free_x", page = page_no)) 
  
}
plot_row = 4
plot_col = 4
pages<- ceiling(cols/(plot_row*plot_col))
for (i in 1:pages){
page_plot(i, plot_row, plot_col)
  
}

2.7 Korelacja miedzy danymi

“Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji.”

# OK
#TODO wiecej wynikow

plot_data <- corr_cross(var_df, 
max_pvalue = 0.05,
top = 100,
rm.na=TRUE)
## Returning only the top 100. You may override with the 'top' argument
#plot(plot_data)
ggplotly(plot_data)

2.8 Zmiana [atrybutu/ow] w czasie

Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie.

timeline_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    layer(data=cleaned, 
          mapping=aes(
              x=re_date, 
              y=hypersensitive_cardiac_troponin_i, group=patient_id
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(timeline_plot)
timeline_plot <- ggplot(cleaned, aes(x=re_date, y=serum_chloride, colour=factor(patient_id), group=patient_id))  + geom_line() + geom_point() + facet_wrap(~outcome)
ggplotly(timeline_plot)
#timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose)
timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose, outcome, gender) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose, outcome, gender) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose), outcome=first(outcome), gender=first(gender)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose, outcome, gender)
## `summarise()` regrouping output by 'id' (override with `.groups` argument)
timeline_plot <- ggplot(timeline_data,aes(day_no,avg_hemoglobin, shape=factor(gender)))+geom_point(color="blue")+geom_point(aes(day_no,avg_glucose, shape=factor(gender)),color="black") + facet_wrap(~outcome)
ggplotly(timeline_plot)
## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.

2.9 Klasyfikator

klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji.

ml_data<- cleaned%>%group_by(patient_id)%>%summarise_all(funs(last))
ml_data<-na_mean(ml_data) #TO ZROBIC OSOBNO DLA ZBIORU TESTUJACEGO I UCZACEGO!
#rozwazyc: usuniecie tych kolumn (wierszy), w których jest dużo wartosci NA (np. powyżej 40%?)

#ml_data$outcome=as.factor(ml_data$outcome)

ml_data$outcome=factor(ml_data$outcome, 
                        labels = make.names(c("negative", "positive")))

set.seed(23)
inTraining <- 
    createDataPartition(
        # atrybut do stratyfikacji
        y = ml_data$outcome,
        # procent w zbiorze uczącym
        p = .75,
        # chcemy indeksy a nie listę
        list = FALSE)
training <- ml_data[ inTraining,]
testing <- ml_data[ -inTraining,]
ctrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5)


fit <- train(outcome ~ .,
             data = training,
             method = "rf",
             trControl = ctrl,
             ntree = 10)
rfClasses <- predict(fit, newdata = testing)
confusionMatrix(data = rfClasses, testing$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction negative positive
##   negative       48        0
##   positive        0       41
##                                                
##                Accuracy : 1                    
##                  95% CI : (0.9594, 1)          
##     No Information Rate : 0.5393               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 1                    
##                                                
##  Mcnemar's Test P-Value : NA                   
##                                                
##             Sensitivity : 1.0000               
##             Specificity : 1.0000               
##          Pos Pred Value : 1.0000               
##          Neg Pred Value : 1.0000               
##              Prevalence : 0.5393               
##          Detection Rate : 0.5393               
##    Detection Prevalence : 0.5393               
##       Balanced Accuracy : 1.0000               
##                                                
##        'Positive' Class : negative             
## 
ml_data$outcome=factor(ml_data$outcome, 
                        labels = make.names(levels(ml_data$outcome)))

rfGrid <- expand.grid(mtry = 10:30)
gridCtrl <- trainControl(
    method = "repeatedcv",
    summaryFunction = twoClassSummary,
    classProbs = TRUE,
    number = 2,
    repeats = 5)

set.seed(23)
fitTune <- train(outcome ~ .,
             data = training,
             method = "rf",
             metric = "ROC",
             preProc = c("center", "scale"),
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 30)
rfTuneClasses <- predict(fitTune,
                         newdata = testing)
confusionMatrix(data = rfTuneClasses, 
                testing$outcome)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction negative positive
##   negative       48        0
##   positive        0       41
##                                                
##                Accuracy : 1                    
##                  95% CI : (0.9594, 1)          
##     No Information Rate : 0.5393               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 1                    
##                                                
##  Mcnemar's Test P-Value : NA                   
##                                                
##             Sensitivity : 1.0000               
##             Specificity : 1.0000               
##          Pos Pred Value : 1.0000               
##          Neg Pred Value : 1.0000               
##              Prevalence : 0.5393               
##          Detection Rate : 0.5393               
##    Detection Prevalence : 0.5393               
##       Balanced Accuracy : 1.0000               
##                                                
##        'Positive' Class : negative             
## 
ggplot(fitTune) + theme_bw()

### Analizę ważności atrybutów najlepszego znalezionego modelu

2.10 Dodatkowa analiza

analiza typowa dla danych klinicznych, np.:

  • regresja logistyczna wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors)
  • regresja Coxa (ang. Cox Proportional-Hazards Model).